
#include <math.h>

#include <util/math_macros.h>
#include <util/msg.h>
#include <util/mem.h>
#include <util/types.h>

#include "_draw_primitives.h"

#define swap( a, b, type ) { type tmp; tmp = a; a = b; b = tmp; }
#define abs( val ) (val)<0?-(val):(val)
#define min( a, b ) ( a < b ? a : b )
#define round(a) ((a)>=0) ?	(long)((a)+0.5)	: (long)((a)-0.5)

static unsigned char * __array;
static int __width;
static int __height;
static int __row_stride;
static int __depth;

static int __color;
void _set_color(int cor) {
	__color = cor;
}
void _set_array(unsigned char * data, int width, int height, int row_stride) {
	__array = data;
	__width = width;
	__height = height;
	__row_stride = row_stride;
	__depth = row_stride / width;
}
inline void _put_pixel(int x, int y) {
	int i = (x + __width * y) * __depth;
	mem_copy((char *) &__array[i], (char *) &__color, __depth);
}
void _bresenham_line(int x0, int y0, int x1, int y1) {

	// = draw vertical line
	if (x0 == x1) {
		if (y0 > y1) swap(y0, y1, int);
		for (; y0 < y1; y0++) _put_pixel(x0, y0);
		return;
	}

	// = draw horizontal line
	if (y0 == y1) {
		if (x0 > x1) swap(x0, x1, int);
		for (; x0 < x1; x0++) _put_pixel(x0, y0);
		return;
	}

	// = draw non vertical and non horizontal line
	int dx = abs(x1 - x0);
	int dy = -abs(y1 - y0);

	int sx = x0 < x1 ? 1 : -1;
	int sy = y0 < y1 ? 1 : -1;

	int err = dx + dy;
	int e2; // error value e_xy

	for (;;) { // loop 

		_put_pixel(x0, y0);

		if (x0 == x1 && y0 == y1) break;
		e2 = 2 * err;
		if (e2 >= dy) {
			err += dy;
			x0 += sx;
		} // e_xy+e_x > 0
		if (e2 <= dx) {
			err += dx;
			y0 += sy;
		} // e_xy+e_y < 0
	}
}
void _bresenham_line_incx(int x0, int y0, int x1, int y1) {

	// = draw vertical line
	if (x0 == x1) {
		if (y0 > y1) swap(y0, y1, int);
		for (; y0 < y1; y0++) _put_pixel(x0, y0);
		return;
	}

	// = draw horizontal line
	if (y0 == y1) {
		if (x0 > x1) swap(x0, x1, int);
		for (; x0 < x1; x0++) _put_pixel(x0, y0);
		return;
	}

	// = draw non vertical and non horizontal line
	int x, y;
	int dx, dy;
	int _2dx, _2dy;
	int d;
	int nSlope;
	int swapped = 0;

	// deal with vertical lines
	if (x0 == x1) {
		if (y0 > y1) swap(y0, y1, int);
		for (y = y0; y <= y1; y++) _put_pixel(x0, y);
		return;
	}
	// deal with horizontal lines
	if (y0 == y1) {
		if (x0 > x1) swap(x0, x1, int);
		for (x = x0; x <= x1; x++) _put_pixel(x, y0);
		return;
	}

	// if dy > dx line is in the 1st octant
	if (abs(y1 - y0) > abs(x1 - x0)) {
		// swap the axes
		swap(x0, y0, int);
		swap(x1, y1, int);
		swapped = 1;
	}

	dx = x1 - x0;
	dy = y1 - y0;

	nSlope = (dy > 0 && dx < 0) || (dy < 0 && dx > 0);
	int sslope = nSlope ? -1 : 1;

	d = sslope * dx - 2 * dy;

	_2dx = sslope * dx << 1;
	_2dy = dy << 1;

	y = y0;
	for (x = x0; x <= x1; x++) {

		if (!swapped) _put_pixel(x, y);
		else _put_pixel(y, x);

		d -= _2dy;

		if (sslope * d < 0) {
			y = y + 1;
			d += _2dx;
		}
	}
}
void _bresenham_line_incx1(int x1, int y1, int x2, int y2) {
	int slope;
	int dx, dy, incE, incNE, d, x, y;
	// Reverse lines where x1 > x2
	if (x1 > x2) {
		_bresenham_line_incx1(x2, y2, x1, y1);
		return;
	}

	dx = x2 - x1;
	dy = y2 - y1;
	// Adjust y-increment for negatively sloped lines
	if (dy < 0) {
		slope = -1;
		dy = -dy;
	} else {
		slope = 1;
	}
	// Bresenham constants
	incE = 2 * dy;
	incNE = 2 * dy - 2 * dx;
	d = 2 * dy - dx;
	y = y1;
	// Blit
	for (x = x1; x <= x2; x++) {

		_put_pixel(x, y);

		//putpixel(x, y);

		if (d <= 0) {
			d += incE;
		} else {
			d += incNE;
			y += slope;
		}
	}
}
void _bresenham_line_incx2(int xl, int yl, int xr, int yr) {
	int x, y; /* coordinates of pixel being drawn  */
	int dy, dx;
	int ne; /* integer scaled error term     */
	int ie;
	x = xl;
	y = yl; /* start at left endpoint            */

	dx = xr - xl;
	dy = yr - yl;

	ie = 2 * dy - dx; /* initialize the error term */
	while (x <= xr) { /* pixel-drawing loop */

		_put_pixel(x, y);
		//PlotPixel(x, y); /* draw the pixel   */
		if (ie > 0) {
			y = y + 1;
			ne = ne - 2 * dx; /* replaces e = e - 1  */
		}
		x = x + 1;
		ne = ne + 2 * dy; /* replaces e = e + m  */
	}
}
void _bresenham_circle(int xm, int ym, int ray) {
	int x = -ray, y = 0, err = 2 - 2 * ray; /* II. Quadrant */
	do {
		_put_pixel(xm - x, ym + y); /*   I. Quadrant */
		_put_pixel(xm - y, ym - x); /*  II. Quadrant */
		_put_pixel(xm + x, ym - y); /* III. Quadrant */
		_put_pixel(xm + y, ym + x); /*  IV. Quadrant */

		ray = err;
		if (ray > x) err += ++x * 2 + 1; /* e_xy+e_x > 0 */
		if (ray <= y) err += ++y * 2 + 1; /* e_xy+e_y < 0 */
	} while (x < 0);
}
void _ellipse_rect(int x0, int y0, int x1, int y1) {
	int a = abs(x1 - x0), b = abs(y1 - y0), b1 = b & 1; /* values of diameter */
	long dx = 4 * (1 - a) * b*b, dy = 4 * (b1 + 1) * a*a; /* error increment */
	long err = dx + dy + b1 * a*a, e2; /* error of 1.step */

	if (x0 > x1) {
		x0 = x1;
		x1 += a;
	} /* if called with swapped points */
	if (y0 > y1) y0 = y1; /* .. exchange them */
	y0 += (b + 1) / 2;
	y1 = y0 - b1; /* starting pixel */
	a *= 8 * a;
	b1 = 8 * b*b;

	do {
		_put_pixel(x1, y0); /*   I. Quadrant */
		_put_pixel(x0, y0); /*  II. Quadrant */
		_put_pixel(x0, y1); /* III. Quadrant */
		_put_pixel(x1, y1); /*  IV. Quadrant */
		e2 = 2 * err;
		if (e2 >= dx) {
			x0++;
			x1--;
			err += dx += b1;
		} /* x step */
		if (e2 <= dy) {
			y0++;
			y1--;
			err += dy += a;
		} /* y step */
	} while (x0 <= x1);

	while (y0 - y1 < b) { /* too early stop of flat ellipses a=1 */
		_put_pixel(x0 - 1, y0); /* -> finish tip of ellipse */
		_put_pixel(x1 + 1, y0++);
		_put_pixel(x0 - 1, y1);
		_put_pixel(x1 + 1, y1--);
	}
}
void _ellipse(int xc, int yc, int a, int b, real deg) {
	float t;
	float inc;
	int x0, x1;
	int y0, y1;

	real cosrad, sinrad;
	real cost, sint;

	if (deg == 0) {
		int a2 = a >> 1;
		int b2 = b >> 1;

		int x0 = xc - a2;
		int x1 = xc + a2;
		int y0 = yc - b2;
		int y1 = yc + b2;

		_ellipse_rect(x0, y0, x1, y1);
		return;
	}

	deg = -(deg * M_PI) / 180;

	inc = 0.01;
	t = 0.0;
	cosrad = cos(deg);
	sinrad = sin(deg);

	cost = cos(t);
	sint = sin(t);

	x0 = xc + a * cost * cosrad - b * sint * sinrad;
	y0 = yc + a * cost * sinrad + b * sint * cosrad;

	t += inc;

	do {
		cost = cos(t);
		sint = sin(t);

		x1 = xc + a * cost * cosrad - b * sint * sinrad;
		y1 = yc + a * cost * sinrad + b * sint * cosrad;

		x0 = x1;
		y0 = y1;
		t += inc;

		_bresenham_line(x0, y0, x1, y1);

	} while (t <= (M_PI + M_PI));

}
void _basic_bezier(int x0, int y0, int x1, int y1, int x2, int y2) {
	int sx = x0 < x2 ? 1 : -1, sy = y0 < y2 ? 1 : -1; /* step direction */
	int cur = sx * sy * ((x0 - x1)*(y2 - y1)-(x2 - x1)*(y0 - y1)); /* curvature */
	int x = x0 - 2 * x1 + x2, y = y0 - 2 * y1 + y2, xy = 2 * x * y * sx*sy;
	/* compute error increments of P0 */
	long dx = (1 - 2 * abs(x0 - x1)) * y * y + abs(y0 - y1) * xy - 2 * cur * abs(y0 - y2);
	long dy = (1 - 2 * abs(y0 - y1)) * x * x + abs(x0 - x1) * xy + 2 * cur * abs(x0 - x2);
	/* compute error increments of P2 */
	long ex = (1 - 2 * abs(x2 - x1)) * y * y + abs(y2 - y1) * xy + 2 * cur * abs(y0 - y2);
	long ey = (1 - 2 * abs(y2 - y1)) * x * x + abs(x2 - x1) * xy - 2 * cur * abs(x0 - x2);

	/* sign of gradient must not change */
	//assert((x0 - x1)*(x2 - x1) <= 0 && (y0 - y1)*(y2 - y1) <= 0);

	if (cur == 0) {
		_bresenham_line(x0, y0, x2, y2);
		return;
	} /* straight line */

	x *= 2 * x;
	y *= 2 * y;
	if (cur < 0) { /* negated curvature */
		x = -x;
		dx = -dx;
		ex = -ex;
		xy = -xy;
		y = -y;
		dy = -dy;
		ey = -ey;
	}
	/* algorithm fails for almost straight line, check error values */
	if (dx >= -y || dy <= -x || ex <= -y || ey >= -x) {
		_bresenham_line(x0, y0, x1, y1); /* simple approximation */
		_bresenham_line(x1, y1, x2, y2);
		return;
	}
	dx -= xy;
	ex = dx + dy;
	dy -= xy; /* error of 1.step */

	for (;;) { /* plot curve */
		_put_pixel(x0, y0);
		ey = 2 * ex - dy; /* save value for test of y step */
		if (2 * ex >= dx) { /* x step */
			if (x0 == x2) break;
			x0 += sx;
			dy -= xy;
			ex += dx += y;
		}
		if (ey <= 0) { /* y step */
			if (y0 == y2) break;
			y0 += sy;
			dx -= xy;
			ex += dy += x;
		}
	}
}
void _bezier(int x0, int y0, int x1, int y1, int x2, int y2) {
	/* plot any quadratic Bezier curve */
	int p, x = x0 - 2 * x1 + x2, y = y0 - 2 * y1 + y2;
	float t = x;
	if ((x0 - x1)*(x2 - x1) > 0) { /* horizontal split at P4? */
		if ((y0 - y1)*(y2 - y1) > 0) /* vertical split at P6 too? */
			if (abs(x0 - x1) * abs(y) > abs(y0 - y1) * abs(x)) { /* which first? */
				p = x0;
				x0 = x2;
				x2 = p; /* if vertical swap points */
				p = y0;
				y0 = y2;
				y2 = p;
			}
		/* now horizontal split at P4 comes first */
		t = (x0 - x1) / t;
		p = (1 - t)*((1 - t) * y0 + 2 * t * y1) + t * t*y2; /* By(t=P4) */
		x = t = (x0 * x2 - x1 * x1) * t / (x0 - x1); /* gradient dP4/dx=0 */
		y = (x0 - t)*(y0 - y1) / (x1 - x0) + y0; /* intersect P3 | P0 P1 */
		_basic_bezier(x0, y0, x, y, x, p);
		y1 = (y1 - y2)*(t - x2) / (x1 - x2) + y2; /* intersect P4 | P1 P2 */
		x0 = x1 = x;
		y0 = p; /* P0 = P4, P1 = P8 */
	}
	if ((y0 - y1)*(y2 - y1) > 0) { /* vertical split at P6? */
		t = y0 - 2 * y1 + y2;
		t = (y0 - y1) / t;
		p = (1 - t)*((1 - t) * x0 + 2 * t * x1) + t * t*x2; /* Bx(t=P6)*/
		y = t = (y0 * y2 - y1 * y1) * t / (y0 - y1); /* gradient dP6/dy=0*/
		x = (t - y0)*(x1 - x0) / (y1 - y0) + x0; /* intersect P6 | P0 P1*/
		_basic_bezier(x0, y0, x, y, p, y);
		x1 = (x1 - x2)*(t - y2) / (y1 - y2) + x2; /* intersect P7 | P1 P2 */
		x0 = p;
		y0 = y1 = y; /* P0 = P6, P1 = P7 */
	}
	_basic_bezier(x0, y0, x1, y1, x2, y2); /* remaining part */

}

typedef struct _edge {
	Ref2 p[2];
	int dx;
	int dy;
	real x;
	real inv_slope;
	int active;

} Edge;
static Ref2 _line_intersec(Ref2 A, Ref2 B, Ref2 C, Ref2 D) {

	Ref2 p;
	real t;

	int a = B.x - A.x;
	int b = D.x - C.x;
	int c = B.y - A.y;
	int d = D.y - C.y;
	int x = C.x - A.x;
	int y = C.y - A.y;

	real u = (real) (a * y - c * x) / (real) (c * b - a * d);

	t = u * b + C.x;
	p.x = ROUND(t);
	t = u * d + C.y;
	p.y = ROUND(t);

	return p;
}
static int _compare_edges(Edge * e1, Edge * e2) {

	int y1, y2;

	y1 = min(e1->p[0].y, e1->p[1].y);
	y2 = min(e2->p[0].y, e2->p[1].y);

	return y1 - y2;
}
static int _compare_x(int * x1, int * x2) {
	return (*x1) - (*x2);
}
void _scan_line_poly_fill(Ref2 * p, int n) {

	int y;
	int ymin, ymax;
	int x_edg_i[4 * n];
	int nic; // number of intersections
	int i;

	Ref2 pi;
	Ref2 p0, p1, p2, p3;

	if (n < 3) return;

	// = Find the upper point and the lower point
	ymin = ymax = p[0].y;
	for (i = 0; i < n; i++) {
		if (p[i].y < ymin) ymin = p[i].y;
		if (p[i].y > ymax) ymax = p[i].y;
	}

	if (ymin == ymax) return;

	for (y = ymin; y <= ymax; y++) {

		// = Search all intersection of horizontal line in y with all edges in the polygon
		nic = 0;

		for (i = 0; i < n; i++) {

			p1 = p[i];
			p2 = p[ (i + 1) == n ? 0 : (i + 1)];

			// = Avoid horizontal edges.
			if (p1.y == p2.y) continue;

			// = Find the intersection with the line  p[i-1] *------* p[i]
			pi = _line_intersec(_Ref2(0, y), _Ref2(1, y), p1, p2);

			// = Verify if the intersections is the vertex p1
			if (pi.x == p1.x && pi.y == p1.y) {

				p0 = p[ (i - 1) < 0 ? (n - 1) : (i - 1)];

				if (p0.y < pi.y && pi.y <= p2.y) {
					// is monotonically decreasing
					if (nic > 0 && x_edg_i[nic - 1] == pi.x) {
						continue;
					}
					x_edg_i[nic] = pi.x;
					nic++;
					continue;
				}
				if (p2.y <= pi.y && pi.y < p0.y) {
					// is monotonically increasing
					if (nic > 0 && x_edg_i[nic - 1] == pi.x) {
						continue;
					}
					x_edg_i[nic] = pi.x;
					nic++;
					continue;
				}

				// is a local minimum or maximum
				if (nic > 0 && x_edg_i[nic - 1] == pi.x) {
					continue;
				}
				x_edg_i[nic] = pi.x;
				nic++;
				x_edg_i[nic] = pi.x;
				nic++;
				continue;

			}
			// = Verify if the intersections is the vertex p2
			if (pi.x == p2.x && pi.y == p2.y) {

				p3 = p[ (i + 2) >= n ? 1 : (i + 2)];

				if (p1.y < pi.y && pi.y <= p3.y) {
					// is monotonically decreasing
					if (nic > 0 && x_edg_i[nic - 1] == pi.x) {
						continue;
					}
					x_edg_i[nic] = pi.x;
					nic++;
					continue;
				}
				if (p1.y <= pi.y && pi.y < p3.y) {
					// is monotonically increasing
					if (nic > 0 && x_edg_i[nic - 1] == pi.x) {
						continue;
					}
					x_edg_i[nic] = pi.x;
					nic++;
					continue;
				}

				// is a local minimum or maximum
				if (nic > 0 && x_edg_i[nic - 1] == pi.x) {
					continue;
				}
				x_edg_i[nic] = pi.x;
				nic++;
				x_edg_i[nic] = pi.x;
				nic++;
				continue;
			}

			//= Verify the intersection with the segment
			if ((p1.y <= pi.y && pi.y <= p2.y) || (p2.y <= pi.y && pi.y <= p1.y)) {
				x_edg_i[nic] = pi.x;
				nic++;
			}
		}
		// = test
		if (nic == 0) continue;
		if (nic % 2 != 0) {
			nic--;
			//msg_warning("Ooops!\n");
		}

		// = Sort all intersections in x direction
		gqsort(x_edg_i, nic, sizeof (x_edg_i[0]), _compare_x);

		for (i = 0; i < nic; i += 2) {
			_bresenham_line(x_edg_i[i], y, x_edg_i[i + 1], y);
		}

	}
}
void _scan_line_poly_fill_fast(Ref2 * p, int n) {
	Edge EL[n];

	int nc;
	int C[n * 4];

	Ref2 p1, p2;
	int i;
	int y;

	int ymin, ymax;
	int ne = 0; // Number of non-horizontal edges

	ymax = p[0].y;

	// = Create all edges and get the ymax
	for (i = 0; i < n; i++) {

		p1 = p[i];
		p2 = p[ (i + 1) == n ? 0 : (i + 1)];

		// = Find ymax
		if (p1.y > ymax) ymax = p1.y;

		// = The first point of each edge has the minimal y coordinate of the two points
		if (p1.y == min(p1.y, p2.y)) {
			EL[i].p[0] = p1;
			EL[i].p[1] = p2;
		} else {
			EL[i].p[0] = p2;
			EL[i].p[1] = p1;
		}

		EL[i].dx = p1.x - p2.x;
		EL[i].dy = p1.y - p2.y;
		EL[i].inv_slope = (real) EL[i].dx / (real) EL[i].dy;
		EL[i].active = 0;

		ne++;
	}

	// = sort edges in relation to the min y coordinate of its points
	gqsort(EL, ne, sizeof (EL[0]), _compare_edges);

	// = Get ymin and ymax coordinates of edges set.
	ymin = EL[0].p[0].y;
	//ymax was calculated above

	int start = 0;

	for (y = ymin; y <= ymax; y++) {

		// = Find the min Edge that intersects with y
		while (!(EL[start].p[0].y <= y && y <= EL[start].p[1].y)) {
			start++;
		}

		nc = 0;
		for (i = start; i < ne; i++) {

			if (y == EL[i].p[0].y) {
				if (y == EL[i].p[1].y) {
					// EL[i] is a horizontal line
					EL[i].active = 0;
					EL[i].x = EL[i].p[1].x;
					C[nc] = EL[i].x;
					nc++;

				} else {
					// EL[i] is a non-horizontal line
					EL[i].active = 1;
					EL[i].x = EL[i].p[0].x;
				}
			}

			if (y == EL[i].p[1].y) {
				// y intersects with the bigger y coordinate of EL[i]
				EL[i].active = 0;
				EL[i].x = EL[i].p[1].x;
				C[nc] = EL[i].x;
				nc++;
			}

			if (EL[i].p[0].y < y && y < EL[i].p[1].y) {
				// y intersects with EL[i]
				EL[i].x += EL[i].inv_slope;
				C[nc] = round(EL[i].x);
				nc++;
			}
		}

		// = sort x coordinates in left to right order
		gqsort(C, nc, sizeof (C[0]), _compare_x);

		if (nc % 2 != 0) {
			msg_warning("Ooops!\n");
		}

		// Draw line
		for (i = 0; i < nc; i += 2) {
			_bresenham_line(C[i], y, C[i + 1], y);
		}

	}
}
int _scan_line_poly_weight_fast(unsigned char * array, int width, Ref2 * p, int n) {

	Edge EL[n];

	int nc;
	int C[n * 4];

	Ref2 p1, p2;
	int i;
	int y;

	int ymin, ymax;
	int ne = 0; // Number of non-horizontal edges

	ymax = p[0].y;

	// = Create all edges and get the ymax
	for (i = 0; i < n; i++) {

		p1 = p[i];
		p2 = p[ (i + 1) == n ? 0 : (i + 1)];

		// = Find ymax
		if (p1.y > ymax) ymax = p1.y;

		// = The first point of each edge has the minimal y coordinate of the two points
		if (p1.y == min(p1.y, p2.y)) {
			EL[i].p[0] = p1;
			EL[i].p[1] = p2;
		} else {
			EL[i].p[0] = p2;
			EL[i].p[1] = p1;
		}

		EL[i].dx = p1.x - p2.x;
		EL[i].dy = p1.y - p2.y;
		EL[i].inv_slope = (real) EL[i].dx / (real) EL[i].dy;
		EL[i].active = 0;

		ne++;
	}

	// = sort edges in relation to the min y coordinate of its points
	gqsort(EL, ne, sizeof (EL[0]), _compare_edges);

	// = Get ymin and ymax coordinates of edges set.
	ymin = EL[0].p[0].y;
	//ymax was calculated above

	int start = 0;
	long int sum = 0;

	for (y = ymin; y <= ymax; y++) {

		// = Find the min Edge that intersects with y
		while (!(EL[start].p[0].y <= y && y <= EL[start].p[1].y)) {
			start++;
		}

		nc = 0;
		for (i = start; i < ne; i++) {

			if (y == EL[i].p[0].y) {
				if (y == EL[i].p[1].y) {
					// EL[i] is a horizontal line
					EL[i].active = 0;
					EL[i].x = EL[i].p[1].x;
					C[nc] = EL[i].x;
					nc++;

				} else {
					// EL[i] is a non-horizontal line
					EL[i].active = 1;
					EL[i].x = EL[i].p[0].x;
				}
			}

			if (y == EL[i].p[1].y) {
				// y intersects with the bigger y coordinate of EL[i]
				EL[i].active = 0;
				EL[i].x = EL[i].p[1].x;
				C[nc] = EL[i].x;
				nc++;
			}

			if (EL[i].p[0].y < y && y < EL[i].p[1].y) {
				// y intersects with EL[i]
				EL[i].x += EL[i].inv_slope;
				C[nc] = ROUND(EL[i].x);
				nc++;
			}
		}

		// = sort x coordinates in left to right order
		gqsort(C, nc, sizeof (C[0]), _compare_x);

		if (nc % 2 != 0) {
			msg_warning("Ooops!\n");
		}

		// Draw line
		int _r, idx;
		int x0;
		int x1;

		for (i = 0; i < nc; i += 2) {

			// = Remove pair copies
			if (C[i] == C[i + 1]) continue;

			x0 = C[i];
			x1 = C[i + 1];
			_r = y * width;

			for (; x0 < x1; x0++) {
				idx = x0 + _r;
				sum += array[idx];
			}
		}
	}
	return sum;
}
int _scan_line_poly_weight_fast2(unsigned char * array, int width, Ref2 * p, int n) {

	Edge EL[n];

	int nc;
	int C[n * 4];

	Ref2 p1, p2;
	int i;
	int y;

	int ymin, ymax;
	int ne = 0; // Number of non-horizontal edges

	ymax = p[0].y;

	// = Create all edges and get the ymax
	for (i = 0; i < n; i++) {

		p1 = p[i];
		p2 = p[ (i + 1) == n ? 0 : (i + 1)];

		// = Find ymax
		if (p1.y > ymax) ymax = p1.y;

		// = The first point of each edge has the minimal y coordinate of the two points
		if (p1.y == min(p1.y, p2.y)) {
			EL[i].p[0] = p1;
			EL[i].p[1] = p2;
		} else {
			EL[i].p[0] = p2;
			EL[i].p[1] = p1;
		}

		EL[i].dx = p1.x - p2.x;
		EL[i].dy = p1.y - p2.y;
		EL[i].inv_slope = (real) EL[i].dx / (real) EL[i].dy;
		EL[i].active = 0;

		ne++;
	}

	// = sort edges in relation to the min y coordinate of its points
	gqsort(EL, ne, sizeof (EL[0]), _compare_edges);

	// = Get ymin and ymax coordinates of edges set.
	ymin = EL[0].p[0].y;
	//ymax was calculated above

	int start = 0;
	long int sum = 0;

	for (y = ymin; y <= ymax; y++) {

		// = Find the min Edge that intersects with y
		while (!(EL[start].p[0].y <= y && y <= EL[start].p[1].y)) {
			start++;
		}

		nc = 0;
		for (i = start; i < ne; i++) {

			if (y == EL[i].p[0].y) {
				if (y == EL[i].p[1].y) {
					// EL[i] is a horizontal line
					EL[i].active = 0;
					EL[i].x = EL[i].p[1].x;
					C[nc] = EL[i].x;
					nc++;

				} else {
					// EL[i] is a non-horizontal line
					EL[i].active = 1;
					EL[i].x = EL[i].p[0].x;
				}
			}

			if (y == EL[i].p[1].y) {
				// y intersects with the bigger y coordinate of EL[i]
				EL[i].active = 0;
				EL[i].x = EL[i].p[1].x;
				C[nc] = EL[i].x;
				nc++;
			}

			if (EL[i].p[0].y < y && y < EL[i].p[1].y) {
				// y intersects with EL[i]
				EL[i].x += EL[i].inv_slope;
				C[nc] = ROUND(EL[i].x);
				nc++;
			}
		}

		// = sort x coordinates in left to right order
		gqsort(C, nc, sizeof (C[0]), _compare_x);

		if (nc % 2 != 0) {
			msg_warning("Ooops!\n");
		}

		// Draw line
		int _r, idx;
		int x0;
		int x1;

		int inc = 3;

		for (i = 0; i < nc; i += 2) {

			// = Remove pair copies
			if (C[i] == C[i + 1]) continue;

			x0 = C[i];
			x1 = C[i + 1];
			_r = y * width;

			idx = x0 + inc + _r;
			sum = array[idx];

			x0 += inc;
			x1 -= inc;
			for (; x0 < x1; x0++) {
				idx = x0 + _r;
				if (sum != array[idx]) return 0;
			}
		}
	}
	return 1;
}